Code for Chapters 1.1 to 1.8
health_full = read_csv("https://chronicdata.cdc.gov/api/views/6vp6-wxuq/rows.csv?accessType=DOWNLOAD")
health_ca = filter(health_full, StateAbbr == "CA")
pge_q1 = read_csv("PGE_2019_Q1_ElectricUsageByZip.csv")
write_csv(health_ca, "health_ca.csv")
saveRDS(health_ca, "health_ca.rds")
save(health_ca, pge_q1, file = "working_datasets.rda")
load("working_datasets.rda")
save.image("progress1.rda")
load("progress1.rda")
year = 2019
quarters = 1:4
type = "Electric"
pge_19_elec = NULL
for(quarter in quarters) {
filename = paste0("PGE_", year, "_Q", quarter, "_", type, "UsageByZip.csv")
print(filename)
temp = read_csv(filename)
pge_19_elec = rbind(pge_19_elec, temp)
saveRDS(pge_19_elec, "pge_19_elec.rds")
}
## [1] "PGE_2019_Q1_ElectricUsageByZip.csv"
## [1] "PGE_2019_Q2_ElectricUsageByZip.csv"
## [1] "PGE_2019_Q3_ElectricUsageByZip.csv"
## [1] "PGE_2019_Q4_ElectricUsageByZip.csv"
pge_filter = filter(pge_19_elec, CUSTOMERCLASS %in% c("Elec- Residential", "Elec- Commercial"))
pge_select = select(pge_filter, -c(YEAR, COMBINED, AVERAGEKWH))
pge_group = group_by(pge_select, MONTH, CUSTOMERCLASS)
pge_summarize = summarize(pge_group, TOTALKWH = sum(TOTALKWH, na.rm = T))
pge_wide = pivot_wider(pge_summarize, names_from = CUSTOMERCLASS, values_from = TOTALKWH)
pge_tidy = pivot_longer(pge_wide, c("Elec- Commercial", "Elec- Residential"), names_to = "CUSTOMERCLASS", values_to = "TOTALKWH")
pge_summarize = summarize(pge_group, TOTALKWH = sum(TOTALKWH, na.rm = T), TOTALCUSTOMERS = sum(TOTALCUSTOMERS, na.rm = T))
pge_mutate = mutate(pge_summarize, AVERAGEKWH = TOTALKWH/TOTALCUSTOMERS)
pge_final = pge_19_elec %>%
filter(CUSTOMERCLASS %in%
c("Elec- Commercial", "Elec- Residential")) %>%
select(!c(YEAR, COMBINED, AVERAGEKWH)) %>%
group_by(MONTH, CUSTOMERCLASS) %>%
summarize(TOTALKWH = sum(TOTALKWH,
na.rm = T),
TOTALCUSTOMERS = sum(TOTALCUSTOMERS,
na.rm = T)) %>%
mutate(AVERAGEKWH = TOTALKWH/TOTALCUSTOMERS)
pge_chart = pge_final %>%
ggplot() +
geom_bar(aes(x = MONTH %>% factor(),
y = TOTALKWH,
fill = CUSTOMERCLASS),
stat = "identity",
position = "stack") +
labs(x = "Month",
y = "kWh",
title = "PG&E Territory Monthly Electricity Usage, 2019",
fill = "Electricity Type")
pge_chart

pge_chart %>% ggplotly() %>%
layout(xaxis = list(fixedrange = T),
yaxis = list(fixedrange = T)) %>%
config(displayModeBar = F)
pge_chart

Code for Chapter 1.9
ca_counties = counties("CA", cb = T, progress_bar = F)
st_crs(ca_counties)
## Coordinate Reference System:
## User input: NAD83
## wkt:
## GEOGCRS["NAD83",
## DATUM["North American Datum 1983",
## ELLIPSOID["GRS 1980",6378137,298.257222101,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["latitude",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["longitude",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4269]]
projection = "+proj=utm +zone=10 +ellps=GRS80 +datum=NAD83 +units=ft +no_defs"
ca_counties_transformed = ca_counties %>%
st_transform(4326) %>%
st_transform(26910) %>%
st_transform(projection) %>%
st_transform(st_crs(ca_counties))
ggplot(ca_counties) + geom_sf()

leaflet() %>% addTiles() %>%
addPolygons(data = ca_counties %>%
st_transform(4326)) %>%
addMarkers(data = ca_counties %>%
st_centroid() %>%
st_transform(4326))
bay_county_names <-
c(
"Alameda",
"Contra Costa",
"Marin",
"Napa",
"San Francisco",
"San Mateo",
"Santa Clara",
"Solano",
"Sonoma"
)
bay_counties = ca_counties %>% filter(NAME %in% bay_county_names)
ggplot(bay_counties) + geom_sf()

ca_cities = places("CA", cb = T, progress_bar = FALSE)
bay_cities = ca_cities[bay_counties, ]
bay_cities_within = ca_cities %>%
st_centroid() %>%
.[bay_counties, ] %>%
st_set_geometry(NULL) %>%
left_join(ca_cities %>% select(GEOID)) %>%
st_as_sf()
leaflet() %>% addTiles() %>%
addPolygons(data = bay_counties %>%
st_transform(4326),
fill = F,
weight = 2) %>%
addPolygons(data = bay_cities %>%
filter(!GEOID %in% bay_cities_within$GEOID) %>%
st_transform(4326),
color = "red") %>%
addPolygons(data = bay_cities_within %>%
st_transform(4326),
color = "green")